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Abstract. We show how to use experiments over finite fields to gain infor- 
mation about the solution set of polynomial equations in characteristic 
zero. 



Introduction 

Let X be a variety defined over Z. According to Grothendieck we can picture X 
as a family of varieties Xp over spec Z with fibers over closed points of spec Z 
corresponding to reductions modulo p and the generic fiber over (0) corresponding 
to the variety Xq defined by the equations of X over Q. 




Co) 

Figure 1. A variety over specZ 



The generic fiber is related to the special fibers by semicontinuity theorems. 
For example, the dimension of Xp is upper semicontinuous with 

dimXo = mindimXo. 

p>0 ^ 

This allows us to gain information about Xq by investigating Xp which is often 
computationally much simpler. 

Even more surprising is the relation between the geometry of Xp and the 
number of Fp rational points of Xp discovered by Weil: 

Theorem 0.1. Let Xp C Pp he a smooth curve of genus g, and N be the number 
of ¥p-rational points of Xp. Then 

\l'N+p\<2g^. 

He conjectured even more precise relations for varieties of arbitrary dimension 
which were proved by Deligne using ^adic cohomology. 

In this tutorial we will use methods which are inspired by Weil's ideas, but are 
not nearly as deep. Rather we will rely on some basic probabilistic estimates which 
are nevertheless quite useful. I have learned these ideas from my advisor Frank 
Schreyer, but similar methods have been used independently by other people, for 
example Joachim von zur Gathen and Igor Shparlinski [1 , Oliver Labs ^ and 
Noam Elkies [5]. 

The structure of these notes is as follows: We start in Section [T] by evaluat- 
ing the polynomials defining a variety X at random points. This can give some 
heuristic information about the codimension c of X and about the number d of 
codimension-c components of X. 

In Section |2] we refine this method by looking at the tangent spaces of X in 
random points. This gives a way to also estimate the number of components in 
every codimension. As an application we show how this can be applied to gain 
new information about the Poincare center problem. 

In Section [3] we explain how it is often possible to prove that a solution found 
over Fp actually lifts to Q. This is applied to the construction of new surfaces in 

Often one would like not only to prove the existence of a lift, but explicitly 
find one. It is explained in Section [4] how this can be done if the solution set is 
zero dimensional. 

We close in Section [5] with a beautiful application of these lifting techniques 
found by Oliver Labs, showing how he constructed a new septic with 99 real nodes 
in Pi 

For all experiments in this tutorial we have used the computer algebra system 
Macaulay 2 [Ij . The most important Macaulay 2 commands used are explained in 
Appendix]^ for more detailed information we refer to the online help of Macaulay 
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2 |3]. In Appendix [b] Stefan Wiedmann provides a MAGMA translation of the 
Macualay 2 scripts in this tutorial. All scripts are available online at [5] . We would 
like to include translations to other computer algebra packages, so if you are for 
example a Singular-expert, please contact us. 

Finally I would like to thank the referee for many valuable suggestions. 



1. Guessing 

We start by considering the most simple case, namely that of a hypersurface 
X C A" defined by a single polynomial / € Fp[xi, . . . ,Xn]- If a G A" is a point 
we have 



/(«) 



one possibility 
7^ (p — 1) possibihties 



Naively we would therefore expect that we obtain zero for about ^ of the points. 

Experiment 1.1. We evaluate a given polynomial in 700 random points, using 
Macaulay 2: 

R = ZZ[x,y,z,w] — work in AA~4 

F = x"23+1248*y*z+w+129269698 — a Polynomial 

K = ZZ/7 — work over F_7 

L = apply (700, — substitute 700 

i->sub(F,rEaid.om(K~l,K~4))) — random points 

tally L — count the results 



obtaining: 



o5 = Tally{-1 => 100} 
-2 => 108 
-3 => 91 

=> 98 

1 => 102 

2 => 101 

3 => 100 



Indeed, all elements of Fy occur about 700/7 = 100 times as one would expect 
naively. 



If f = g ■ h Cz ¥,p[xi, . . . , Xn] is a reducible polynomial we have 
/(a) = g{a)h{a) 



0-01 possibility 
* • (p — 1) possibilities 
• * (p — 1) possibilities 



(p— 1) possibilities 



so one might expect a zero for about ^^2^ ~ | "-"^ the points. 

Experiment 1.2. We continue Experiment |1.1| and evaluate a product of two 
polynomials in 700 random points: 

G = x*y*z*w+z~25-938493+x-z*w — a second polynomial 
tally apply (700, — substitute 700 

i->sub(F*G, random (K~1,K~4))) — random points & count 



This gives: 



08 = Tally-C-l => 86} 
-2 => 87 
-3 => 77 

=> 198 

1 => 69 

2 => 84 

3 => 99 



Indeed, the value now occurs about twice as often, i.e. 198 ~ | • 700. 

Repeating Experiments |1.1| and |1.2| for 100 random polynomials and 100 
random products we obtain Figure [2j 
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Figure 2. Evaluating 100 random polynomials and 100 random products at 700 points each. 

Observe that the results for irreducible and reducible polynomials do not 
overlap. Evaluating a polynomial at random points might therefore give some 
indication on the number of its irreducible factors. For this we will make the above 
naive observations more precise. 
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Definition 1.3. If / G Fp[a;i, . . . , x„] is a polynomial, we call the map 

/If," : F^' ^ F, 
a /(a) 

the corresponding polynomial function. We denote by 

Vp:={f:¥;^¥,} 
the vector space of all polynomial functions on F^. 

Being a polynomial function is nothing special: 

Lemma 1.4 (Interpolation). Let (f>: Fp — > Fp be any function. Then there exists a 
polynomial f € Vp[xi, . . . , Xn] such that (f> — f\fri . 

Proof. Notice that (1 — a;^^^) = x ^ 0. For every a £ F^ we define 

n 

fa{x) ■.= l[{l-{x,-a,r-') 
1=1 

and obtain 

ft \ _/ lifa; = a 
Ja[x) - |oif x^a. 

Since F^ is finite we can consider / :— J2ae¥" 't'{^)fa ^^^d obtain f{x) = (t>{x) for 
all x e F^. " □ 

Remark 1.5. From Lemma Fl .41 it follows that 

(i) is a vector space of dimension p". 

(ii) is a finite set with pf elements. 

(iii) Two distinct polynomials can define the same polynomial function, for ex- 
ample x^ and X. More generally if _F: Fp — s- Fp is the Frohenius endomor- 
phism then /(a) = f{F{a)) for all polynomials / and all a G F^. 

This makes it easy to count polynomial functions: 
Proposition 1.6. The number of polynomial functions / G with k zeros is 

^'J^ .lfc.(p_l)p"-^ 

Proof. Since Vp is simply the set of all functions /: F^ — > Fp, we can enumerate 
the ones with k zeros as follows: First choose k points and assign the value and 
then chose any of the other (p — 1) values for the remaining — k points. □ 
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Corollary 1.7. The average number of zeros for polynomial functions f Vp is 
and the standard deviation of the number of zeros in this set is 



"fG)(^)<^- 

Proof. Standard facts about binomial distributions. □ 

Remark 1.8. Using the normal approximation of the binomial distribution, we 
can estimate that more than 99% of all / e satisfy 

- /i| < 2.58^^ 

For products of polynomials we have 

Proposition 1.9. The number of pairs {f,g)(zVpXVp whose product has k zeros 
is 

• .9) - fc} = i^l) ■ (2p- 1)'^ • {{p-lfY'^-K 
In particular, the average number of zeros in this set is 

and the standard deviation is 



Proof. As in the proof of Proposition |1.6| we first choose k points. For each of 
these points x we choose either the value of f{x) = and g{x) ^ or f{x) ^ 
and g(x) = or f{x) = g{x) = 0. This gives 2p— 1 possibilities. For the remaining 
p" — k we choose / and q nonzero. For this we have {p — 1)^ possibilities. The 
formulas then follow again from standard facts about binomial distributions. □ 

Remark 1.10. It follows that more than 99% of pairs {f,g)^VpXVp satisfy 

\W{f-9)-l^'\<2.58^. 

In particular, if a polynomial / has a number of zeros that lies outside of this 
range one can reject the hypothesis that / is a product of two irreducible with 
99% confidence. 
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Figure 3. Distribution of the number of zeros on hypersurfaces in A'* in characteristic 7. 



Even for small p the distributions of Proposition [L6] and Proposition |1.9| difFer 
substantially (see Figure |3]). 

Remark 1.11. For plane curves we can compare our result to the Weil conjectures. 
Weil shows that 100% of smooth pane curves of genus g in P^^ satisfy 

\N-{p+l)\<2g^ 

while we proved that 99% of the polynomial functions on satisfy 

|iV-p| <2.58V^. 

Of course Weil's theorem is much stronger. If p > Ag"^ Weil's theorem implies for 
example that every smooth curve of genus g over Fp has a rational point, while 
no such statement can be derived from our results. If on the other hand one is 
satisfied with approximate results, our estimates have the advantage that they 
are independent of the genus g. In Figure |4] we compare the two results with an 
experiment in the case of plane quartics. (Notice that smooth plane quartics have 
genus 3.) 

For big n it is very time consuming to count all Fp-rational points on V{f) C 
A". We can avoid this problem by using a statistical approach once again. 

Definition 1.12. Let X C A" be a variety over Fp. Then 

#^(Fp) 



#A«(Fp) 
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Figure 4. Number of points on 1000 alEne quartics in characteristic g = 37 compared to the 
corresponding binomial distribution and Weil's bound. 

is called the fraction of zeros of X. If furthermore ii, . . . , € A"(Fp) are points 
then 

7p X 

m 

is called an empirical fraction of zeros of X. 

Remark 1.13. If we choose the points Xi randomly and independently, the prob- 
ability that Xi e X is jp{X). Therefore we have the following: 

(i) /i(7p) = 7p, i.e. for large m we expect %{X) w JpiX) 

(ii) cr(7p) « w i.e. the quadratic mean of the error |7p — jp\ decreases with 



/m. 

(iii) Since for hypersurfaces 7p « ^ the average error depends neither on the 
number of variables n nor on the degree of X. 

(iv) Using the normal approximation again one can show that it is usually 
enough to test about 100 • p points to distinguish between reducible and 
irreducible polynomials (for more precise estimates see [B]). 

Experiment 1.14. Consider quadrics in P'^ and let 

A := {singular quadric} C {all quadrics} = 

be the subvariety of singular quadrics in the space of all quadrics. Since having 
a singularity is a codimension 1 condition for surfaces in P"^ we expect A to be a 
hypersurface. Is A irreducible? Using our methods we obtain a heuristic answer 
using Macaulay 2: 
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— work in characteristic 7 
K = ZZ/7 

— the coordinate Ring of IP~3 
R = K[x,y,z,w] 

— look at 700 quadrics 

tally apply(700, i->codim singularLocus (ideal rajidom(2,R) ) ) 

giving 

ol2 = Tally{2 => 5 } 

3 => 89 

4 => 606. 

We see 95 = 89 + 5 of our 700 quadrics were singular, i.e. 7(A) — Since this 
is much closer to then it is to | we guess that A is irreducible. Notice that we 
have not even used the equation of A to obtain this estimate. 

Let's now consider an irreducible variety X C A" of codimension c > 1. 
Projecting A" to a subspace A"^^+^ we obtain a projection X' C A"""^"*"^ of 
X (see Figure |5|. Generically X' is a hypersurface, so by our arguments above 
X' has approximately p"""^ points. Generically most points of X' have only one 
preimage in X so we obtain the following very rough heuristic: 

Heuristic 1.15. Let X C AJ be a variety of codimension c and d the number of 
components of codimension c, then 



Remark 1.16. A more precise argument for this heuristic comes from the Weil 
Conjectures. Indeed, the number of Fp-rational points on an absolutely irreducible 
projective variety X is 

pdimX _|_ j^Q^g]. order terms, 

so 7p « pcodim X ■ Our elementary arguments still work in the case of complete 
intersections and determinantal varieties [3]. 

Remark 1.17. Notice that Heuristic 11.151 involves two unknowns: c and d. To 
determine these one has to measure over several primes of good reduction. 



Experiment 1.18. As in Experiment 1.14 we look at quadrics in P''. These are 



given by their 10 coefficients and form a P^. This time we are interested in the 
variety X C P^ of quadrics whose singular locus is at least one dimensional. For 
this we first define a function that looks at random quadrics over Fp until it has 
found at least k examples whose singular locus has codimension at most c. It then 
returns the number of trials needed to do this. 
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/A 




Figure 5. The projection of a curve in is a hypersurface in and most points have only one 
preimage. 



findk = (p,k,c) -> ( 
K := ZZ/p; 
R : = K [x , y , z , w] ; 
trials := 0; 
f ound : = ; 
while found < k do ( 
Q := ideal raiidoin(2,R) ; 
if c>=codim (Q+ideal jacobian Q) then ( 
found = foimd + 1; 
print found; 

); 

trials = trials + 1; 

); 

trials 
) 



Here we use (Q+ideal jacobian Q) instead of singularLocus(Q), since the 
second option quickly produces a memory overflow. 

The function findk is useful since the error in estimating 7 from 7 depends 
on the number of singular quadrics found. By searching until a given number of 
singular quadrics is found make sure that the error estimates will be small enough. 

We now look for quadrics that have singularities of dimension at least one 
k=50; time LI = apply({5,7, ll},q->(q,time f indk(q,k,2))) 
obtaining 
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Figure 6. Measuring tlie codimension of quadrics with zero and one dimensional singular loci. 
Here we compare the measurements with lines of the correct slope 1 and 3. 

{(5, 5724), (7, 17825), (11, 68349)} 

i.e. 75 ~ sfii; 77 ~ Vf§25 ^'^^ ~ 6^49' '^^^'^ codimension c of X can be 
interpreted as the negative slope in a log-log plot of 'jp{X) since Heuristic 1.15 
gives 

%{X) ~ 4 ^ log(7p(^)) « log('^) - clog(p). 

This is illustrated in Figure [6] 

By using f indk with fc = 50 the errors of all our measurements are of the 
same magnitude. We can therefore use regression to calculate the slope of a line 
fitting these measurements: 

— calculate slope of regression line by 

— formula from [2] p. 800 
slope = (L) -> ( 

xbar := sum (apply (L,1->1#0))/#L; 
ybar := sum(apply(L,l->l#l))/#L; 
sum(apply (L , l-> (l#0-xbar) * (l#l-ybar) ) ) / 
sum(apply(L,l->(l#0-xbar) ~2)) 

) 

— slope for dim 1 singularities 
slope(apply(Ll,l->(log(l/l#0) ,log(k/l#l)))) 

o5 = 3.13578 
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The codimension of X is indeed 3 as can be seen by the following geometric 
argument: Each quadric with a singular locus of dimension 1 is a union of two 
hyperplanes. Since the family P-^ of all hyperplanes in P'^ is 3-dimensional, we 
obtain dimX = 6 which has codimension 3 in the P^ of all quadrics. 

The approach presented in this section measures the number of components 
of minimal codimension quite well. At the same time it is very difficult to see 
components of larger codimension. One reason is that the rough approximations 
that we have made introduce errors in the order of —hrr. 

We will see in the next section how one can circumvent these problems. 



2. Using Tangent Spaces 

If X C A" has components of different dimensions, the guessing method of Section 
[ijdoes not detect the smaller components. 

If for example X is the union of a curve and a surface in A^, we expect the 
surface to have about points while the curve will have about p points (see 
Figure [7| . 




Figure 7. Expected number of ¥p rational points on a union of a curve and a surface. 



Using Heuristic 1.15 we obtain 

_ p2 + p ^ 1 
P"^ p 

indicating that X has 1 component of codimension 1. The codimension 2 compo- 
nent remains invisible. 

Experiment 2.1. Let's check the above reasoning in an experiment. First define a 
function that produces a random inhomogeneous polynomial of given degree: 

rEmdomAf f ine = (d,R) -> sum apply(d+l,i->random(i,R)) 
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with this we choose random polynomials F, G and in 6 variables 



n=6 

R=ZZ[x_l . .x_n] ; 
F = randomAff ine(2,R) 
G = rsindomAf f ine (6 ,R) ; 
H = rcUidomAf f ine (7 ,R) ; 

and consider the ideal / = {FG, FH) 

I = ideal (F*G,F*H) ; 

Finally, we evaluate the polynomials of / in 700 points of characteristic 7 and 
count how many of them lie in X = V{I): 

K = ZZ/7 

t = tally apply (700, i->( 

== sub(I,random(K~l,K'ii)) 
)) 

This yields 

o9 = Tally{false => 598} 
true => 102 

i.e. "fT^X) — which is very close to i. Consequently we would conclude that 
X has one component of codimension 1. The codimension 2 component given by 
G = H = remains invisible. 



To improve this situation we will look at tangent spaces. Let a e X C A" be 
a point and Tx,a the tangent space of X in a. If Ix = {fi, ■ ■ ■ , fm), let 



/ d;^^ 



Jx = 



dfm 



dfi 

dXn 



dXn 



be the Jacobian matrix. We know from differential geometry that 

Tx,a = ker Jx{a) = e if " | Jx{a)v = 0}. 
We can use tangent spaces to estimate the dimension of components of X: 

Proposition 2.2. Let a Cz X d A" be a point and X' <Z X a component containing 
a. Then dimX' < dimTx.a with equality holding in smooth points of X . 



Proof. ^ II. 1.4. Theorem 3] 



□ 
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In particular, we can use the dimension of the tangent space in a point a £ X 
to separate points that he on different dimensional components, at least if these 
components are non reduced (see Figure [s]). For each of these sets we use Heuristic 
II. 151 to obtain 




Figure 8. Dimension of tangent spaces in Fp rational points on a union of a curve and a surface. 

Heuristic 2.3. Let X C A" be a variety . If Jx is the Jacobian matrix of X and 
ai, . . . , am € A" are points, then the number of codimension c components of X 
is approximately 

#{i \ai G X and rank Jxi^i) — c} ■ p"^ 
m 



Experiment 2.4. Let's test this heuristic by continuing Experiment 2.1 For this 
we first calculate the Jacobian matrix of the ideal / 



J = jacobicin I; 



Now we check again 700 random points, but when we find a point on X — V{I) 
we also calculate the rank of the Jacobian matrix in this point: 



K=ZZ/7 

time t = tally apply (700, i->( 
point := random(K~l ,K~n) ; 
if subd, point) == then 
rank sub(J, point) 

)) 



The result is 



ol2 = Tally{0 => 2 > 

1 => 106 

2 => 14 
null => 578 
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Indeed, we find that there are about "'"^^'J — 1.06 components of dimension 1 and 
about ^^'J^ =0.98 components of codimension 2. For codimension the result is 
~ 0.003 consistent with the fact that there are no components of codimension 

0. 



Remark 2.5. It is a little dangerous to give the measurements as in Experiment |2.4| 
without error bounds. Using the Poisson approximation of binomial distributions 
with small success probability we obtain 



tj(number of points found) « i/number of points found. 

In the above experiment this gives 

,. , (106 ±2.58^106) , „^ 

# codim 1 components = = 1.06 ± 0.27 

700 

and 

f 14 ± 2 58\/l4l • 7^ 

# codim 2 components = ^ '- — - — = 0.98 ± 0.68. 

^ ^ 700 

where the error terms denote the 99% confidence interval. Notice that the mea- 
surement of the codimension 2 components is less precise. As a rule of thumb 
good error bounds are obtained if one searches until about 50 to 100 points of 
interest are found. 

Remark 2.6. This heuristic assumes that the components do not intersect. If com- 
ponents do have high dimensional intersections, the heuristic might give too few 
components, since intersection points are singular and have lower codimensional 
tangent spaces. 

In more involved examples calculating and storing the Jacobian matrix Jx 
can use a lot of time and space. Fortunately one can calculate Jx{o,) directly 
without calculating Jx first: 

Proposition 2.7. Let f G ¥p[xi, . . . , Xn] be a polynomial, a G ¥p a point and 
b € ¥p a vector. Then 

f{a + be) = f{a) + 4/(a)e e F^[£]/(£2). 

with djjf denoting the derivative of f in direction ofb. In particular, if ei G is 
the i-th unit vector, we have 

f{a + e^e) = /(a) + ^(a)^- 

Proof. Use the Taylor expansion. □ 
Example 2.8. fix) = =^ /(I + e) ^ (l + e)^ ^ I + 2e ^ /(I) -I- £/'(l) 
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Experiment 2.9. To compare the two methods of calculating derivatives, we con- 
sider the determinant of a random matrix with polynomial entries. First we create 
a random matrix 



calculate the determinant 

time F = det M; 

— used 13.3 seconds 

and its derivative with respect to xi. 

time Fl = diff(x_l,F); 

— used 0.01 seconds 

Now we substitute a random point: 

point = random(K~l ,K"6) 

time sub(Fl .point) 

— used 0. seconds 

o7 = 2 

By far the most time is used to calculate the determinant. With the £-method 
this can be avoided. We start by creating a vector in the direction of xi: 

T = K[e]/(e'2) — a ring with e"2=0 

el = matrix{-[l, 0,0, 0,0,0}} — the first unit vector 

pointl = sub(point,T) + e*sub(el,T) — point with direction 

Now we first evaluate the matrix M in this vector 

time Ml = sub(M, pointl) 

— used 0. seconds 

and only then take the determinant 

time det sub (M, pointl) 

— used 0. seconds 

ol2 = 2e + 1 



K = ZZ/7 

R = K[x_l. .x_6] 

M = random(R~{5:0},R"{5:-2}) 



characteristic 7 

6 variables 

a random 5x5 matrix with 
quadratic entries 
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Figure 9. A focus and a center. 

Indeed, the coefficient of e is the derivative of the determinant in this point. 

This method is too fast to measure by the time command of Macaulay 2. To get 
a better time estimate, we calculate the derivative of the determinant at 5000 
random points: 

time apply (5000, i->( 

point := rcaidom(K~l ,K~6) ; — random point 

pointl := sub(point,T)+e*sub(el,T) ; — tangent direction 
det sub (M, point 1) ; — calculate derivative 

)); 

— used 12.76 seconds 

Notice that this is still faster than calculating the complete determinant once. 



Remark 2.10. The ^-method is most useful if there exists a fast algorithm for 
evaluating the polynomials of interest. The determinant of an n x n matrix for 
example has n\ terms, so the time to (n'aluate it directly is proportional to n!. If 
we use Gauss elimination on the matrix first, the time needed drops to rv^. 

For the remainder of this section we will look at an application of these 
methods to the Poincare center problem. We start by considering the well known 
system of differential equations 

x= -y 
y = x 

whose integral cnirves are circles around the origin. Let's now disturb these equa- 
tions with polynomials P and Q whose terms have degree at least 2: 

X = -y + P 
y = x + Q. 
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Figure 10. Geometric interpretation of the components of tlie center variety in the case d = 2. 

Near zero the integral curves of the disturbed system are either closed or not. In 
the second case one says that the equations have a focus in (0, 0) while in the first 
case they have a center (see Figure [9]). 

The condition of having a center is closed in the space of all (P, Q): 

Theorem 2.11 (Poincare). There exists an infinite series of polynomials fi in the 
coefficients of P and Q such that 



x = -y + 
y = X + ' 



has a center 



f,iP,Q) = Ofor alii. 



We call fi{P, Q) the i-th focal value of (P, Q). 

If the terms of P and Q have degree at most d then the fi describe an algebraic 
variety Xao in the finite- dimensional space of pairs {P,Q). This variety is called 
the center variety. 

Remark 2.12. By Hilbert's Basis Theorem loo {fo, fi, ■ ■ ■) is finitely gener- 
ated. Unfortunately, Hilbert's Basis Theorem is not constructive, so it is a priory 
unknown how many generators lac has. It is therefore useful to consider the i-th 
partial center varieties Xi = V{fo, . . . , fi). 

The following is known: 
Theorem 2.13. If d — 2 then the center variety has four components 

X^=XhU Xiii U Xij UXiC A^, 

three of codimension 2 and one of codimension 3. Moreover X^c = X^. 



Proof. Decompose I3 = (/i, /2, /a) with a computer algebra system and show that 
all solutions do have a center [5], [TU]. □ 
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• p=17, 7 equations, 80 billion points 

• p=19, 8 equations, 90 billion points 

• p=23, 10 equations, 330 billion points 
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Figure 11. Measurements for the Poincare center problem with d - 



Looking at algebraic integral curves one even obtains a geometric interpreta- 



tion of the components in this case (see Figure 10) 



For d~'i almost nothing is known. The best results so far are lists of centers 
given by Zoladec[llj. [T^j . The problem from a computer algebra perspective is 
that the /i are too large to be handled, already /s has 5348 terms and it is known 
that ^ Xio. 

Experiment 2.14. Fortunately for our method, Frommer [9] has devised an algo- 
rithm to calculate /i (P, Q) for given (P, Q) . A closer inspection shows that From- 
mer's Algorithm works over finite fields and will also calculate fi{P+eP' , Q+eQ'). 



So we have all ingredients to use Heuristic 2.3 Using a fast C-|— I- implementation 
of Frommer's Algorithm by Martin Cremer and Jacob Kroker [13j we first check 
our method on the known degree 2 case. For this we evaluate /i, . . . , /lo for d = 2 
at 1.000.000 random points in characteristic 23. This gives 
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88 



Heuristic 12.31 translates this into 



codim components: 0.00 +/- 0.00 
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Table 1. Known families of cubic centers that could have codimension below 8 in A^* |12 |. 



Type 


Name 


Codimension 


Darboux 


CDi 


5 


Darboux 


CD2 


6 


Darboux 


CDs 


7 


Darboux 


CD4 


7 


Darboux 


CDs 


7 


Reversible 


CRi 


> 6 


Reversible 


CR5 


> 7 


Reversible 


CRt 


> 7 


Reversible 


CRu 


> 7 


Reversible 


CR12 


> 7 


Reversible 


CR16 


> 7 



codim 1 components: 0.00 +/- 0.00 
codim 2 components: 2.87 +/- 0.10 
codim 3 components: 1.07 +/- 0.29 



This agrees well with Theorem |2.11[ 

For d = 3 we obtain the measurements in Figure 11 One can check these 
results against Zoladec's lists as depicted in Table[l] Here the measurements agree 
in codimension 5 and 6. In codimension 7 there seem to be 8 known families while 
we only measure 4. Closer inspection of the known families reveals that CRq and 
CRj are contained in CD4 and that CR12 and CRie are contained in CD2 [14] . 
After accounting for this our measurement agrees with Zoladec's results and we 
conjecture that Zoladec's lists are complete up to codimension 7. 



3. Existence of a Lift to Characteristic Zero 



Often one is not interested in characteristic p solutions, but in solutions over C. 
Unfortunately, not all solutions over Fp lift to characteristic 0. 

Example 3.1. Consider the variety X ~ V{3x) C over specZ. As depicted in 
Figure 12 X decomposes into two components: V{3) = Pj.^ which lives only over 
F3 and V (x) = {(0 : 1)} which has fibers over all of specZ. In particular, the 
point (1 : 0) e Pj.^ C X does not lift to characteristic 0. 

To prove that a given solution point over Fp does lift to characteristic zero 
the following tool is very helpful: 

Proposition 3.2 (Existence of a Lifting). Let X C Y C be varieties with 
dimYf^ = dhnYz — 1 for all p and X CY determinantal, i.e. there exists a vector 
bundle morphism 

6: E-^ F 
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Figure 12. The vanishing set of 3x in Pg over spec Z 

on Y and a number r < min(ranki?,ranki^) such that X = Xr{4>) is the locus 
where (j) has rank at most r. If x € Xf^ is a point with 

dim Txpp ,x — dim Ip^ — (rank E — r) (rank F — r) 

then X is smooth in x and there exists a component Z of X^ containing x and 
having a nonzero fiber over (0). 

Proof. Set d = dimYpp — (rankii' — r)(ranki^ — r). Since X^^ is determinantal, 
we have 

dim Zf^ > d 

for every irreducible component Zf of Xf and d is the expected dimension of 
[13 Ex. 10.9, p. 245]. If Zf^ contains the point x we obtain 

d < dim < dim Tz^^ < dim Txj^ = d 

by our assumptions. So Z^^ is of dimension d and smooth in x. Let now Z^ be a 
component of X^ that contains Z^^ and x. Since Xz is determinantal in Yz and 
dim Yz = dim Y^^ + 1 we have 

dim Zz > d+1. 

Since dimZp^ = d the fiber of Z-^ over p cannot contain all of Z-^. Indeed, in this 
case we would have Zf^ — Zi since both are irreducible, but dim Z^^ ^ dim Zz- It 
follows that Zi has nonempty fibers over an open subset of spec Z^ and therefore 
also over (0) [TB], [17]. □ 
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Figure 13. Tangent spaces in several points oi X = V{3x) C over specZ 

Example 3.3. The variety X — V{3x) is determinantal on Y = ¥^ since it is the 
rank locus of the vector bundle morphism 

(f>: Opi —> C'pi(l). 



1 = dim — 1 for all p. The expected dimension of Xr 



- 0) • (1 — 0) = 0. As depicted in Figure 13 we have three typical 



Furthermore dim 
is therefore 1 — (1 
examples: 



(i) a; = (0 : 1) over Fp with p 3. Here the tangent space over ¥p is zero 
dimensional and the point lifts according to Proposition 3.2 

(ii) X = {0 : 1) over F3. Here the tangent space is 1-dimensional and Proposi- 
tion [3^ does not apply. Even though the point does lift. 

(iii) X = {1 : 0) over F3. Here the tangent space is also 1-dimensional and 
Proposition |3 . 2 1 does not apply. In this case the point does not lift. 



This method has been used first by Frank Schreyer |16j to construct new 
surfaces in P** which are not of general type. The study of such surfaces started in 
1989 when EUingsrud and Peskine showed that their degree is bounded [TH] and 
therefore only finitely many families exist. Since then the degree bound has been 
sharpened by various authors, most recently by [l9] to 52. On the other hand a 
classification is only known up to degree 10 and examples are known up to degree 
15 (see |19j for an overview and references). 

Here I will explain how Cord Erdenberger, Katharina Ludwig and I found a 
new family of rational surfaces S of degree 11 and sectional genus 11 in P'' with 
finite field experiments. 
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Our plan is to realize S* as a blowup of P^. First we consider some restrictions 
on the linear system that embeds S into P**: 

Proposition 3.4. Let S = P'^{pi, . . . ,pi) be the blowup o/P^ in I distinct points. 
We denote by Ei,...,Ei the corresponding exceptional divisors and by L the 
pullback of a general line in P^ to S. Let \aL — ^i^il ^ very ample linear 
system of dimension four and set 13 j — #{i | bi — j}. Then 



d^a'-J2f3jf 

j 




= 9-^/3,. 

j 

where d is the degree, tt the sectional genus and K the canonical divisor of S. 

Proof. Intersection theory on S [XT', Corollary 4.1]. □ 

By the double point formula for surfaces in P^ [201 Appendix A, Example 
4.1.3] a rational surface of degree 11 and sectional genus 11 must satisfy — 
— 11. For fixed a the equations above can be solved by integer programming, using 
for example the algorithm described in Chapter 8 of [2T| . 

In the case a < 9 we find that there are no solutions. For a = 9 the only 
solution is /33 = 1, /32 = 14 and /3i = 5. Our first goal is therefore to find 5 simple 
points, 14 double points and one triple point in P^ such that the ideal of the union 
of these points contains 5 polynomials of degree 9. 

To make the search fast, we would like to use characteristic 2. The difficulty 
here is that P^ contains only 7 rational points, while we need 20. Our solution to 
this problem was to choose 

PeP^(F2) Qep2(F2i4) i?GP^(F25) 

such that the Frobenius orbit of Q and R are of length 14 and 5 respectively. The 
ideals of the orbits are then defined over F2 . 

— define coordinate ring of P~2 over F_2 
F2 = GF(2) 

S2 = F2[x,y,z] 

— define coordinate ring of P~2 over F_2~14 and F_2~5 
St = F2[x,y,z,t] 

use St; 114 = ideal(t'14+t'13+t~ll+t~10+t~8+t'6+t"4+t+l) ; S14 = St/114 
use St; 15 = ideal(t"5+t"3+t~2+t+l) ; S5 = St/15 
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— the random points 

use S2; P = matrix{{0_S2 , 0_S2, 1_S2}} 

use S14;Q = matrix{{t'' (random(2"14-l) ) , t~ (random(2~14-l) ) , 1_S14}} 
use S5; R = matrix-[{t" (random 31), t*(random 31), 1_S5}} 

— their ideals 

IP = ideal ((vars S2)*syz P) 

IQ = ideal ((vars S14) _{0 . . 2}*syz Q) 

IR = ideal ((vars S5)_{0. .2}*syz R) 

— their orbits 

fl4 = map (314/10, S2) ; Qorbit = ker fl4 
degree Qorbit — hopefully degree = 14 

f5 = map(S5/IR,S2) ; Rorbit = ker f5 
degree Rorbit — hopefully degree = 5 

If Q and R have the correct orbit length we calculate \QH — 3P —2Q — R\ 

— ideal of 3P 
P3 = IP"3; 

— orbit of 20 

fl4square = map(S14/IQ"2,S2) ; Q2orbit = ker fl4square; 

— ideal of 3P + 2C!orbit + IRorbit 
I = inter sect (P3,Q2orbit, Rorbit ) ; 

— extract 9-tics 

H = super basis (9, I) 

rank source H — hopefully affine dimension = 5 

If at this point we find 5 sections, we check that there are no unassigned base 
points 

— count basepoints (with multiplicities) 

degree ideal H — hopefully degree = 1x6+14x3+1x5 = 53 

If this is the case, the next difficulty is to check if the corresponding lineax system 
is very ample. On the one hand this is an open condition, so it should be satisfied 
by most examples, on the other hand we are in characteristic 2, so exceptional loci 
can have very many points. An irreducible divisor for example already contains 
approximately half of the rational points. 

— construct map to P'"4 
T = F2[x0,xl,x2,x3,x4] 
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fH = map(S2,T,H) ; 



— calculate the ideal of the image 
Isurface = ker fH; 

— check invar icoits 

betti res coker gens Isurface 
codim Isurface — codim = 2 
degree Isurface — degree = 11 
genera Isurface — genera = {0,11,10}- 

— check smoothness 

J = jacobian Isurface; 

mJ = minors (2, J) + Isurface; 

codim mJ — hopefully codim = 5 



Indeed, after about 100.000 trials one comes up with the points 



use S14;q = matrix{{t~ 11898, t~137, 1_S14}} 
use S5; R = matrix{{t~6, t"15, 1_S5» 



These satisfy all of the above conditions and prove that rational surfaces of degree 
11 and sectional genus 11 in exist in over F2. 

As a last step we have to show that this example lifts to char 0. For this we 
consider the morphism 

Tk : H\Or.{a)) Op2(a) © 30p2(a - 1) © • • • © (^^ + ^^Op2(a - fc) 

on P| that associates to each polynomial of degree a the coefficients of its Taylor 
expansion up to degree A; in a given point P. 

Lemma 3.5. If a > k then the image of Tk is a vector bundle J-'k of rank C^^^) 
over specZ. 



Proof. In each point we consider an affinc 2-dimensionaI neighborhood where we 
can choose the (''J^) coefficients of the affine Taylor expansion independently. 
This shows that the image has at least this rank everywhere. If follows from the 
Euler relation for homogeneous polynomials 



df df df 

""di+yTy+'d-z^^'^'^f^-^ 

that this is also the maximal rank. □ 
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Now set = Hilbi^z x Hilbi4,z x Hilbs^z where Hilb^^z denotes the Hilbert 
scheme of k points in P| over specZ, and let 

= {{p, q, r) I h"{9L ~3p~2q- Ir) > 5} C Fz 

be the subset where the hnear system of nine-tics with a triple point in p, double 
points in q and single base points in r is at least of projective dimension 4. 

Proposition 3.6. There exist vector bundles E and F of ranks 55 and 53 respec- 
tively on Yz and a morphism 

(j): E-> F 

such that X5(){(f)) — X^- 

Proof. On the Cartesian product 

Hilbd,zxP| ^ P| 

TTl 

Y 

Hilbd.z 

we have the morphisms 

Let now C Hilb^.z xP^ be the universal set of points. Then P^ is a flat family 
of degree d over Hilb^^z and 

Gk (^i)*((7r;^fe)|pj 
is a vector bundle of rank dC'^^) over Hilb^ z- On 

Fz = Hilbi,z X Hilbi4^z x Hilbs^z 

the induced map 

0: i/°(Op|(9)) ® Ox, Ii^Zl^ ^*g^ ^ a*^g^ ^ ^5*^0 

has the desired properties, where ad denotes the projection to Hilb^.z- □ 

So we have to show that the tangent space of in our base locus has 
codimension (55 — 50) (53 — 50) = 15. This can be done by explicitly calculating 
the differential of cf) in our given base scheme using the e-method. The script is 
too long for this paper, but can be downloaded at (25]. Indeed, we find that the 
codimension of the tangent space is 15, so this shows that our example lies on an 
irreducible component that is defined over an open subset of specZ. 
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speciality 



Figure 14. The difficulty of finding a surface grows exponentially with the speciality 

Remark 3.7. The overall time to find smooth surfaces that lift to characteristic 
zero can be substantially reduced if one calculates the tangent space of a given 
point {P,Q,R) in the Hilbert scheme Xz directly after establishing \9H — 3P — 
2Q — R\ = P^. One then needs to check very ampleness only for smooth points 
of Xz. This is useful since the tangent space calculation is just a linear question, 
while the check for very ampleness requires Grobner bases. We use a very fast C- 
implementation by Jakob Kroker to do the whole search algorithm up to checking 
smoothness. Only the (very few) remaining examples are then checked for very 
ampleness using Macaulay 2. 



Experiment 3.8. We also tried to reconstruct the other known rational surfaces 
in P"' with our program. The number of trials needed is depicted in Figure |14[ 
The expected codimension of X in the corresponding Hilbert scheme turns out to 
be 5 times the speciality h^{Ox{^)) of the surface. As expected, the logarithm of 
the number of trials needed to find a surface is proportional to the codimension 
ofX. 

Remark 3.9. We could not reconstruct all known families. The reason for this 
is that we only look at examples where the base points of a given multiplicity 
form an irreducible Frobenius orbit. In some cases such examples do not exist for 
geometric reasons. 

Experiment 3.10. Looking at the linear system |14iJ — 4P — SQ — i?| with degP = 
8, degQ = 6 and degi? = 2, we find rational surfaces of degree 12 and sectional 
genus 12 in P'* with this method (not published) . 



4. Finding a Lift 

In some good cases characteristic p methods even allow one to find a solution over 
Q quickly. Basically this happens when the solution set is zero dimensional with 
two different flavors. 
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The first good situation, depicted in Figure |15[ arises when X — 
V{fi, . . . , fm) C has a unique solution over Q, maybe with high multiphcity. 
In this case it follows that the solution is defined over O. 




Figure 15. A scheme over spec Z with a unique solution over Q, possibly with high multiplicity 



Algorithm 4.1. If the coordinates of the unique solution over Q are even in Z one 
can find this solution as follows: 

(i) Reduce mod pi and test all points in F^. 

(ii) Find many primes pi with a unique solution in F^. 

(iii) Use Chinese remaindering to find a solution mod JliPi >> 0. 

(iv) Test if this is a solution over Z. If not, find more primes pi with unique 
solutions over Fp . . 

Remark 4.2. Even if the solution y over Q is unique, there can be several solutions 
over Fp. Since the codimension of points in A" is n we expect that the probability 
of a random point x € A" to satisfy x e X is 4r by Heuristic 



expect that the probability oi x ^ X for all points x ^ y is 



1.15 



We therefore 



1 



p" 



Experiment 4.3. Let's use this Algorithm 4.1 to solve 

-Sx^ -xy- 7y2 + 5238x - 11582y - 7696 = 
Axy - lOy^ - 2313x - 16372?; - 6462 = 

For this we need a function that looks at all points over a given prime: 
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allPoints = (I,p) -> ( 
K = ZZ/p; 

flatten apply (p,i-> 
flatten apply (p,j-> 

if (O==codim sub(I,inatrix{{i*l_K, j*l_K}})) 
then {(i.j)} 
else {} 

)) 

) 

With this we look for solutions of our equations over the first nine primes. 

R = ZZ[x,y] 

— the equations 

I = ideal (-8*x"2-x*y-7*y'2+5238*x-11582*y-7696, 
4*x*y-10*y~2-2313*x-16372*y-6462) 

— look for solutions 

tally apply ({2 ,3,5,7,11,13,17,19, 23} , p-> (p , t ime allPoints (I , p) ) ) 

We obtain: 

o8 = Tally{(2, {(0, 0)}) => 1 

(3, {(0, 2), (1, 0), (2, 0)» => 1 

(5, {(4, 1)» => 1 

(7, {(2, 3), (5, 5)}) => 1 

(11, {(2, 7), (8, 1)}) => 1 

(13, {(3, 4), (12, 6)» => 1 

(17, {(10, 8)}) => 1 

(19, {(1, 3), (1, 17), (18, 5), (18, 18)}) => 1 
(23, {(15, 8)}) => 1 

As expected for the intersection of two quadrics we find at most 4 solutions. Over 
four primes we find unique solutions, which is reasonably close to the expected 
number 9/e « 3.31. We now combine the information over these four primes using 
the Chinese remainder Theorem. 

— Chinese remaindering 

— given solutions mod m cind n find 

— a solution mod m*n 

— soli = (n, solution) 

— sol2 = (m, solution) 
chinesePair = (soli, sol2) -> ( 

n = soll#0;an = soll#l; 
m = sol2#0;ain = sol2#l; 
drs = gcdCoeff icients(n,m) ; 



29 



— returns {d,r,s} so that a*r + b*s is the 

— greatest common divisor d of a and b. 
r = drs#l; 

s = drs#2; 

emin = s*m*£m.+r*n*am; 

Eunn = amn - (round(£min/(m*n) ) ) * (m*n) ; 

if (drs#0) == 1 then (m*n,ajiin) else print "m and n not coprime" 
) 

— take a list { (n_l , s_l) , . . . , (n_k, s_k)} 

— and return (n,a) such that 

— n=n_l* ... * n_k and 

— s_i = a mod n_i 

chineseList = (L) -> (fold(L,chinesePair)) 

— X coordinate 

chineseList({(2,0) , (5,4) , (17,10) , (23,15)}) 

— y coordinate 

chineseList({(2,0) , (5,1) , (17,8) , (23,8)}) 

This gives 

011 = (3910, 1234) 

012 = (3910, -774) 

i.e (1234, —774) is the unique solution mod 3910 = 2 • 5 • 17 • 23. Substituting this 
into the original equations over Z shows that this is indeed a solution over Z. 

sub (I , matrix{{1234 , -774}}) 

013 = ideal (0, 0) 

If the unique solution does not have Z but Q coordinates then one can find 
the solution using the extended Euclidean Algorithm [23J Section 5.10]. 

Excimple 4.4. Let's try to find a small solution to the equation 

- = 7 mod 37. 

s 

Each solution satisfies 

r = 7s + 37t 

with s and t in Z. Using the extended Euclidean Algorithm 
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r 


s 


t 


r/s 




37 





1 




-5 


7 


1 





7/1 


-3 


2 


-5 


1 


-2/5 




1 


16 


-3 


1/16 



we find the solution 

1 = gcd(7, 37) = 7 • 16 + 37 • (-3) 

to our linear equation. Observe, however, that the intermediate step in the Eu- 
clidean Algorithm also gives solutions, most of them with small coefficients. In- 
deed, r/s = —2/5 is a solution with r,s < -x/SY which is the best that we can 
expect. 



If we find a small solution by this method, we even can be sure that it is the 
only one satisfying the congruence: 

Proposition 4.5. There exist at most two solutions (r, s) of 

r = as + bt mod m 
that satisfy r,s < \/rn. If a solution satisfies r,s < -ky/m, then this solution is 



2 

unique. 

Proof f53i, Section 5.10] □ 
Experiment 4.6. Let's find a solution to 

176x^ + U8xy + 301y^ - 742a; 4- 896?/ 768 = 
-25xy + 430y^ + 33x + 1373?/ -t- 645 = 



As in Experiment 4.3 we search for primes with unique solutions 



I = ideal (176*x'2+148*x*y+301*y~2-742*x+896*y+768, 
-25*x*y+430*y"2+33*x+1373*y+645) 

tally apply ({2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41}, 
p->(p,time allPoints (I ,p) ) ) 



and obtain 



olO = Tally{(2, {(1, 0)» => 1 

(3, {(0, 0), (0, 1), (2, 0)}) => 1 

(5, {(3, 2), (4, 1)» => 1 

(7, {(2, 6), (4, 0)» => 1 

(11, {» => 1 

(13, {(5, 10)» => 1 
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(17, {(5, 4), (9, 13), (11, 16), (12, 12)}) => 1 

(19, {(3, 15), (8, 6), (13, 15), (17, 1)» => 1 

(23, {(15, 18), (19, 12)}) => 1 

(29, {(26, 15), (28, 9)}) => 1 

(31, {(7, 22)}) => 1 

(37, {(14, 18)}) => 1 

(41, {(0, 23)}) => 1 

Notice that there is no solution mod 11. If there is a solution over Q this means 
that 11 has to divide at least one of the denominators. Chinese remaindering gives 
a solution mod 2 • 13 • 31 • 37 • 41 = 1222702: 

— X coordinate 

chineseList({(2,l) , (13,5) , (31,7) , (37,14) , (41,0)}) 

011 = (1222702, 138949) 

— y coordinate 

chineseList({(2,0) , (13,10) , (31,22) , (37,18) , (41,23)}) 

012 = (1222702, -526048) 

Substituting this into the original equations gives 

sub (I ,matrix{{ 138949 , -526048}}) 

013 = ideal (75874213835186, 120819022681578) 

so this is not a solution over 1. To find a small possible solution over Q we use 
an implementation of the extended Euclidean Algorithm from ^23j Section 5.10]. 

— take (a,n) and calculate a solution to 

r = as mod n 

— such that r,s < sqrt(n). 

— return (r/s) 
recoverQQ = (a,n) -> ( 

rO:=a;sO:=l;tO:=0; 

rl:=n;sl:=0;tl:=l; 

r2:=0;s2:=0;t2:=0; 

k := round sqrt (rl*l . 0) ; 

while k <= rl do ( 

q = rO//rl; 

r2 = rO-q*rl; 

s2 = sO-q*sl; 

t2 = tO-q*tl; 

— print (q,r2,s2,t2) ; 

rO=rl;sO=sl;tO=tl; 

rl=r2;sl=s2;tl=t2; 

); 

(r2/s2) 
) 
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This yields 



— X coordinate 

recoverQQ (138949, 2*13*31*37*41) 
123 

021 = 

22 

Notice that Macaulay reduced 246/44 to 123/22 in this case. Therefore this is not 
a solution mod 2. Indeed, no solution mod 2 exists, since the denominator of the 
X coordinate is divisible by 2. For the y coordinate we obtain 

— y coordinate 

recoverQQ (-526048 , 2*13*31*37*41) 
77 

022 = 

43 

As a last step we substitute this Q-point into the original equations. 

sub (I, matrix{{123/22,-77/43>}) 
o24 = ideal (0, 0) 

This shows that we have indeed found a solution over Q. Notice also that as 
argued above one of the denominators is divisible by 2 and the other by 11. 

Remark 4.7. 

(i) The assumption that we have a unique solution over Q is not as restrictive 
as it might seem. If we have for example 2 solutions, then at least the line 
through them is unique. More generally, if the solution set over Q lies on k 
polynomials of degree d then the corresponding point in the Grassmannian 
G(fc, ('^+")) is unique. 

(ii) Even if wc do not have isolated solutions, we can use this method to find the 
polynomials of rad(/(X)), at least if the polynomials are of small degree. 

(iii) For this method we do not need explicit equations, rather an algorithm that 
decides whether a point lies on X is enough. This is indeed an important 
distinction. It is for example easy to check whether a given hypersurface is 
singular, but very difficult to give an explicit discriminant polynomial in 
the coefficients of the hypersurface that vanishes if and only if it is singular. 

Before we finish this tutorial by looking at a very nice application of this 
method by Oliver Labs, we will look briefly at a second situation in whic;h we can 
find explicit solutions over Q. I learned this method from Noam Elkies in his talk 
at the Clay Mathematics Institute Summer School "Arithmetic geometry" 2006. 
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Figure 16. A scheme X over spec Z with a smooth isolated solution over Q. 

Algorithm 4.8. Assume that X has a smooth point x over Q that is isolated 
over Q as depicted in Figure 16 and that p is a prime that does not divide the 
denominators of the coordinates of x. Then we can find this point as follows: 

(i) Reduce mod p and test all points. 

(ii) Calculate the tangent spaces at the found points. If the dimension of such 
a tangent space is then the corresponding point is smooth and isolated. 

(iii) Lift the point mod p'' with k large using p-adic Newton iteration, as ex- 



plained in Prop 4.9 



Proposition 4.9. Let a E be a solution of 

/i(a)---- = /n(a) = mod/ 
and assume that the Jacobian matrix J — ^^^^ is invertible at a mod p. Then 

a' ^a - (/i(a), . . . , /„(a)) J(a)"^ 

is a solution mod p^^ . 

Proof. Use the Taylor expansion as in the proof of Newton iteration. 



□ 



Experiment 4.10. Let's solve the equations of Experiment 4.3 using p-adic Newton 
iteration. For this we need some functions for modular calculations: 

— calculate reduction of a matrix M mod n 
modn = (M,n) -> ( 

matrix apply (reink target M, i-> 

applyCraoik source M,j-> M_j_i-round(M_j_i/n) *n) ) ) 



— divide a matrix of integers by eoi integer 
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— (in our application this division will not have a remainder) 
divn = (M,n) -> ( 

matrix apply (rank target M, i-> 

apply (rcink source M,j-> M_j_i//n))) 

— invert number mod n 
invn = (i,n) -> ( 

c := gcdCoef f icients (i ,n) ; 

if c#0 == 1 then c#l else "error" 

) 

— invert a matrix mod n 

— Ma square matrix over ZZ 

— (if M is not invertible mod n, then is returned) 
invMatn = (M,n) -> ( 

Mn := modn(M,n) ; 
MQQ := sub(Mn,QQ) ; 
detM = sub(det Mn.QQ); 

modn(invn(sub(detM,ZZ) ,n)*sub(detM*MQQ~-l ,ZZ) ,n) 
) 

With this we can implement Newton iteration. We will represent a point by a 
pair (P, eps) with P a matrix of integers that is a solution modulo eps. 

— (P.eps) an approximation mod eps (contains integers) 

— M affine polynomials (over ZZ) 

— J Jacobian matrix (over ZZ) 

— returns an approximation (P,eps~2) 
newtonStep = (Peps, M, J) -> ( 

P := Peps#0; 
eps := Peps#l; 

JPinv := invMatn(sub(J,P) ,eps) ; 

correction := eps*modn(divn(sub(M,P)*JPinv,eps) ,eps) ; 

■Cmodn(P-correction,eps~2) , eps"2} 

) 

— returns an approximation mod Peps~(2"num) 
newton = (Peps.M, J.num) -> ( 

i := 0; 

localPeps := Peps; 
while i < num do ( 

localPeps = newtonStep(localPeps,M, J) ; 

print (localPeps) ; 

i = i+1; 

); 

localPeps 



35 



We now consider equations of Example 



4.3 



I = ideal (-8*x~2-x*y-7*y"2+5238*x-11582*y-7696 , 
4*x*y-10*y~2-2313*x-16372*y-6462) 



their Jacobian matrix 

J = jacobicind) 

and their solutions over F7: 

apply(allPoiiits(I,7) ,Pseq -> ( 
P := matrix {toList Pseq}; 
(P,0!=det modn(sub(J,P) ,7)) 

)) 

o25 = {(I 2 3 I, true), (| 5 5 I, true)} 



Both points are isolated and smooth over F7 so we can apply p-adic Newton 



iteration to them. The first one lifts to the solution found in Experiment 4.3 



newton((matrix{{2,3}},7) ,gens I, J, 4) 

{| 9 10 I , 49} 

{I -1167 -774 I , 2401} 

{| 1234 -774 I , 5764801} 

{I 1234 -774 I, 33232930569601} 

while the second point probably does not lift to Z: 

newton((matrix{{5,5}},7) ,gens I, J, 4) 

{| 5-9 I , 49} 

{| -926 334 I , 2401} 

{| 359224 -66894 I, 5764801} 

{I 11082657337694 -9795607574104 I, 33232930569601} 

Remark 4.11. Noam Elkies has used this method to find interesting elliptic fibra- 
tions over Q. See for example [31 Section III, p. 11]. 

Remark 4.12. The Newton method is much faster than lifting by Chinese remain- 
dering, since we only need to find one smooth point in one characteristic. Unfor- 
tunately, it does not work if we cannot calculate tangent spaces. An application 
where this happens is discussed in the next section. 
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Figure 17. Historic plaster model of the Cayley Cubic as displayed in the mathematical instute 
of the university of Gottingen 



5. Surfaces with Many Real Nodes 



A very nice application of finite field experiments with beautiful characteristic 
zero results was done by Oliver Labs in his thesis [2]- We look at his ideas and 
results in this section. 

Consider an algebraic surface X C Pj^ of degree d and denote by N{X) the 
number of real nodes oi X. A classical question of real algebraic geometry is to 
determine the maximal number of nodes a surface of degree d can have. We denote 
this number by 



fi{d) := ma.x{N{X) \ X cVl/\degX ^ d}. 



Moreover one would like to find explicit equations for surfaces X that do have 
IJ,{d) real nodes. The cases = and /i(2) = 1, i.e the plane and the quadric 
cone, have been known since antiquity. 

Cayley [24j and Schafli [2S] solved /x(3) = 4, while Kummer proved fi{A) = 16 
in [23] . Plaster models of a Cayley-Cubic and a Kummer-Quartic are on display 
in the Gottingen Mathematical Institute as numbers 124 and 136, see Figure [T7| 



and 18 These pictures many other are available at 



http: / /www.uni-math.gwdg.de/modellsammlung. 



For the case d = 5, Togliatti proved in [27 that quintic surfaces with 31 nodes 
exist. One such surface is depicted in Figure [20] It took 40 years before Beauville 
[28] finally proved that 31 is indeed the maximal possible number. 

In 1994 Barth [59] found the beautiful sextic with the icosahedral symmetry 



and 65 nodes shown in Figure 21 Jaffe and Rubermann proved in [30^ that no 



sextics with 66 or more nodes exist. 
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For d = 7 the problem is still open. By works of Chmutov [3J, 
Breske/Labs/van Straten [32 and Varichenko [33] we only know 93 < /i(7) < 104. 
For large d Chmutov and Breske/Labs/van Straten show 

5 

fi{d) — Y2^'^ ~^ lower order terms, 
while Miyaoka [33] proves 

fJ-id) < gii^ + lower order terms. 

Here we explain how Oliver Labs found a new septic with many nodes, using 
finite field experiments |35j . 

Experiment 5.1. The most naive approach to find septics with many nodes is to 
look at random surfaces of degree 7 in some small characteristic: 

— Calculate milnor number for hypersurf aces in IP~3 

— (for nonisolated singularities sind. smooth surfaces is returned) 
mu = (f) -> ( 

J := (ideal jacobian ideal f)+ideal f; 
if 3==codim J then degree J else 
) 

K = ZZ/5 — work in char 5 

R = K[x,y,z,w] — coordinate ring of IP~3 

— look at 100 random surfaces 

time tally apply (100, i-> mu (random (7 ,R) ) ) 

After about 18 seconds we find 

o4 = Tally{0 => 69> 

1 => 24 

2 => 5 

3 => 1 

4 => 1 

which is still far from 93 nodes. Since having an extra node is a codimension-one 
condition, a rough estimation gives that we would have to search 5®^ « 1.6 x 10^^ 
times longer to find 89 more nodes in characteristic 5. 

One classical idea to find surfaces with many nodes, is to use symmetry. If 
for example we only look at mirror symmetric surfaces, we obtain singularities in 
pairs, as depicted in Figure [19] 
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Figure 18. A Kummor surface with 16 nodes. 

Experiment 5.2. We look at 100 random surfaces that are symmetric with respect 
to the X ~ plane 

— make a random f mirror symmetric 
sym = (f) -> f +sub(f ,{x=>-x}) 

time tally applyClOO, i-> mu(sym(rEaidom(7,R)))) 

o6 = Tally-CO => 57} 

1 => 10 

2 => 11 

3 => 9 

4 => 4 

5 => 3 

6 => 3 

7 => 1 
9 => 1 
13 => 1 

Indeed, we obtain more singularities, but not nearly enough. 

The symmetry approach works best if we have a large symmetry group. In 
the d — 7 case Oliver Labs used the symmetry of the 7-gon. If Dj acts on 
P'^ with symmetry axis x — y = one can use representation theory to find a 
7-dimensional family of _D7-invariant 7-tic we use in the next experiment. 

Experiment 5.3. Start by considering the cone over a 7-gon given by 




which can be expanded to 
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Figure 19. A mirror symmetric cubic 



P = x*(x~6-3*7*x'4*y"2+5*7*x~2*y~4-7*y~6)+ 

7*z*((x~2+y'2)"3-2~3*z~2*(x~2+y~2)"2+2~4*z~4*(x~2+y~2))- 
2-6*z-7 

Now parameterize Di invariant septics U that contain a double cubic. 

S = K[al,a2,a3,a4,a5,a6,a7] 

RS = R**S — tensor product of rings 

U = (z+a5*w)* 

(al*z~3+a2*z~2*w+a3*z*w~2+a4*w"3+(a6*z+a7*w) * (x~2+y~2) ) ~2 

We will look at random sums of the form P + U using 

randomlnv =()->( 

P-sub(U,vars R I random (R~{0>,R"{7:0>)) 
) 

Let's try 100 of these 

time tally apply (100, i-> muCreindomlnvO ) ) 

o9 = Tally{63 => 48} 

64 => 6 

65 => 4 

136 => 1 
140 => 1 

Unfortunately, this looks better than it is, since many of the surfaces with high 
Milnor numbers have singularities that arc not ordinary nodes. We can detect this 
by looking at the Hessian matrix which has rank > 3 only at smooth points and 
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Figure 20. A Togliatti quintic 



ordinary nodes. The following function returns the number of nodes of X = V{f) 
if all nodes are ordinary and otherwise. 



numAl = (f) -> ( 

— singularities of f 

singf := (ideal jacobian ideal f)+ideal f; 
if 3==codim singf then ( 
— calculate Hessian 

Hess := dif f (trsmspose vars R,diff(vars R,f)); 

ssf := singf + minors (3 , Hess) ; 

if 4==codiin ssf then degree singf else 

) 

else 



With this we test another 100 examples: 



time tally applyClOO, i-> numAl(randomInv() ) ) 

ol2 = Tally{0 => 28 > 

63 => 51 

64 => 13 

65 => 1 
70 => 6 
72 => 1 



which takes about 30 seconds. Notice that most surfaces have N{X) a multiple 
of 7 as expected from the symmetry. 



To speed up these calculations Oliver Labs intersects the surfaces X = V{P+ 
U) with the hyperplane y = see Figure 22 Since the operation of moves 
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Figure 21. The Barth sextic 

this hyperplane to 7 different positions, every singularity of the intersection curve 
C that does not he on the symmetry axis corresponds to 7 singularities of X. 
Singular points on C that do lie on the symmetry axis contribute only one node 
to the singularities of X. Using the symmetry of the construction one can show 
that for surfaces X with only ordinary double points all singularities are obtained 
this way [Ml P- 18, Cor. 2.3.10], [Ml Lemma 1]. 

Experiment 5.4. We now look at 10000 random i)7-invariant surfaces and their 
intersection curves with y = Q. We estimate the number of nodes on X from the 
number of nodes on C and return the point in the parameter space of U if this 
number is large enough. 

use R 

time tally apply (10000 , i-> ( 
r := random (R~{0},R" {7:0}) ; 
f := sub(P-sub(U,vars R|r) ,y_R=>0) ; 
singf := ideal f + ideal jacobian ideal f; 
if 2 == codim singf then ( 
— calculate Hessian 

Hess := dif f (trsinspose vars R,diff(vars R,f)); 
ssf := singf + minors (2 , Hess) ; 
if 3==codim ssf then ( 

d := degree singf; 

— points on the line x=0 

singf X := singf +ideal(x) ; 

dx := degree singf x; 

if 2!=codim singf x then dx=0; 

d3 = (d-dx)*7+dx; 

(d,d-dx,dx,d3,if d3>=93 then r) 

) 

else -1 

) 
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Figure 22. Intersection of the 7-gon with a perpendicular hypersurface 

)) 

In this way we find 

ol6 = Tally{(9, 9, 0, 63, ) => 5228 
(10, 9, 1, 64, ) => 731 



(16, 14, 2, 100, I 1221101 I)=>1 
-1 => 3071 
null => 8 

It remains to check whether the found U really gives rise to surfaces with 100 
nodes 

f = P-sub(U,vars R I sub(matrix{-[l ,2, 2 , 1 , 1 , , 1>} ,R) ) 
numAKf ) 

ol8 = 100 

This proves that there exists a surface with 100 nodes over F5. 

Looking at other fields one finds that F5 is a special case. In general one only 
finds surfaces with 99 nodes. To lift these examples to characteristic zero, Oliver 
Labs analyzed the geometry of the intersection curves of the 99-nodal examples 
and found that 

(i) All such intersection curves decompose into a line and a 6-tic. 

(ii) The singularities of the intersection curves are in a special position that 
can be expUcitly described (see [3S] for details) 
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Figure 23. The Labs septic 



These geometric properties imply (after some elimination) that there exists an a 
such that 

ai = + 7a^ - + 7a^ -2o? -la-\ 
a2 = {a^ + l)(3a^ + Ua^ - ia^ + 7a - 2) 
ag = + l)2(3a3 + 7q,_3) 

a" 

"5 = 



1 + a2 
aa — a-j — 1 

It remained to determine which a lead to 99-nodal septics. Experiments over many 
primes show that there are at most 3 such a. Over primes with exactly 3 solutions, 
Oliver Labs represented them as zeros of a degree 3 polynomial. By using the 
Chinese remaindering method, he lifted the coefficients of this polynomial to 
characteristic and obtained 

Ta^ + 7a + 1 = 0. 

This polynomial has exactly one real solution, and with this a one can calculate 



this time over Q(a) that the resulting septic has indeed 99 real nodes. Figure 23 
shows the inner part of this surface. 

A movie of this and many other surfaces in this section can by found on my 
home page 



www.iag. uni-hannover . de / " bothmer /goettingen. php 



on the home page of Oliver Labs 
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http: / /www.algebraicsurface.net / 



or on youTube.com 



http://www.youtube.com/profile?user=botlimer 



The movies and the surfaces in this article were produced using the pubhc domain 
programs surf by Stefan Endrafi [37] and surf ex by Oliver Labs |38l . 



A. Selected Macaulay Commands 

Here we review some Macaulay 2 commands used in this tutorial. Lines starting 
with "i" are input lines, while lines starting with "o" are output lines. For more 
detailed explanations we refer to the online help of Macaulay2 [4^. 

A.l. apply 

This command applies a function to a list. In Macaulay 2 this is often used to 
generate loops. 

11 : apply({l,2,3,4},i->i-2) 
ol = {1, 4, 9, 16} 

01 : List 

The list {0,l,...,n — 1} can be abbreviated by n: 

12 : apply(4,l->l~2) 

02 = {0, 1, 4, 9} 
o2 : List 

A . 2. map 

With map(R,S,m) a map from S" to i? is produced. The matrix m over S contains 
the images of the variables of R: 

11 : f = map(ZZ,ZZ[x,y] ,matrlx{{2,3») 
ol = map(ZZ,ZZ[x,y] ,{2, 3}) 

01 : RlngMap ZZ < ZZ[x,y] 

12 : f(x+y) 

02 = 5 

If no matrix is given, all variables to variables of the same name or to zero. 
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13 : g = map(ZZ [x] ,ZZ [x,y] ) 

o3 = map(ZZ[x] ,ZZ[x,y] ,{x, 0}) 

03 : RingMap ZZ[x] < ZZ[x,y] 

14 : g(x+y+l) 

04 = x + 1 
o4 : ZZ[x] 

A. 3. random 

This command can be used either to construct random matrices 

11 : K = ZZ/3 
ol = K 

01 : QuotientRing 

12 : randoiii(K~2,K"3) 

02 = I 1 -II 

11-11 I 

2 3 

02 : Matrix K < K 

or to construct random homogeneous polynomials of given degree 

13 : R = K[x,y] 

03 = R 

03 : PolynomialRing 

14 : randoiii(2,R) 

2 2 

04 = X + x*y - y 
o4 : R 

A. 4- sub 

This command is used to substitute values for the variables of a ring: 
il : K = ZZ/3 
ol = K 

ol : QuotientRing 
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12 : R = K[x,y] 
o2 = R 

02 : PolynomialRing 

13 : f = x*y 

03 = x*y 

03 : R 

14 : sub(f ,matrix{-[2,3}}) 

04 = 6 

Another application is the transfer a polynomial, ideal or matrix from one ring 
R to another ring 5 that has some variables in common with R 

15 : S = K[x,y,z] 

05 = S 

05 : PolynomialRing 

16 : sub(f,S) 

06 = x*y 
06 : S 

A. 5. syz 

The command is used here to calculate a presentation for the kernel of a matrix: 

11 : M = matrix{{l,2,3},-C4,5,6» 

ol = I 12 3 1 
14 5 6 1 

2 3 

01 : Matrix ZZ < ZZ 

12 : syz M 

02 = I -1 I 

I 2 I 
I -1 I 

3 1 
o2 : Matrix ZZ < ZZ 
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A. 6. tally 



With tally one can count how often an element appears in a hst: 



il : tally[l,2,l,3,2,2,17} 



ol = Tally-Cl => 2 } 

2 => 3 

3 => 1 
17 => 1 

ol : Tally 



B. Magma Scripts (by Stefan Wiedmann) 

Stefan Wiedmann [5^ has translated the Macaulay 2 scripts of this article to 
Magma. Here they are: 

Experiment B.1.1. Evaluate a given polynomial in 700 random points. 



K := FiniteField(7) ; 



//work over F_7 



R<x,y,z,w> := PolynomialRing(K,4) ; //Polynomialring in 4 variables over F_7 



K4:=CartesianPower(K,4) ; 

F := x~23+1248*y*z*w+129269698; 



//K-4 

//a polynomial 



M := [Random (K4) : i in [1..700]] 
T := {*Evaluate(F,s) : s in M*}; 



//random points 



Multiplicity(T,0) ; 



//Results with muliplicity 



Experiment B.1.2. Evaluate a product of two polynomials in 700 random points 



K := FiniteField(7) ; 

R<x,y,z,w> := PolynomialRing(K,4) : 

K4:=CartesianPower(K,4) ; 



//work over F_7 
//AA~4 over F_7 
//K-4 



= x~23+1248*y*z*w+129269698; 
= x*y*z*w+z~25-938493+x-z*w; 
= F*G; 



//a polynomial 

//a second polynomial 



M := [Random (K4) : i in [1..700]]; //random points 

T := {*Evaluate(H,s) : s in M*}; 

T; 

Multiplicity(T,0) ; 

Experiment B.1.14. Count singular quadrics. 



48 



K := FiniteField(7) ; 

R<X , Y , Z , W> : = PolynomialRing (K , 4) ; 

{* Dimension(JacobianIdeal(Randoin(2,R,0))) : i in [1..700]*}; 
Experiment B.1.18. Count quadrics with dim > singular locus 

function f indk(n,p,k,c) 

//Search until k singular examples of codim at most c are found, 

//p prime number, n dimension 

K := FiniteField(p) ; 

R := PolynomialRing (K, n) ; 

trials := 0; 

foimd := 0; 

while foimd It k do 

Q := IdeaK [Random (2,R,0)] ) ; 

if c ge n - Dimension (Q+Jacobicoildeal (Basis (Q))) then 

foimd := found + 1; 
else 

trials := trials + 1; 
end if; 
end while; 

print "Trails: " ,trials; 
return trials; 
end function; 

k := 50; 

time LI := [[p,f indk(4,p,k,2)] : p in [5,7,11]]; 
LI; 

time findk(4,5,50,2) ; 
time findk(4,7,50,2) ; 
time findk(4, 11,50,2) ; 

function slope (L) 

//calculate slope of regression line by 

//formula form [2] p. 800 

xbar := &+[L[i] [1] : i in [l..#L]]/#L; 

ybar := &+[L[i][2] : i in [l..#L]]/#L; 

return &+[(L[i] [l]-xbar)*(L[i] [2]-ybar) : i in [1..3]]/ 
&+[(L[i] [l]-xbar)-2 : i in [1..3]]; 
end function; 

//slope for dim 1 singularities 

slope ([ [Log (l/x[l]), Log(k/x[2])] : x in LI] ) ; 

Experiment B.2.1. Count points on a reducible variety. 
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K := FiniteField(7) ; 

V := CartesicaiPower(K,6) ; 

R<xl,x2,x3,x4,x5,x6> := PolynomialRing(K,6) ; 

//random affine polynomial of degree d 

randomAf f ine := func< d I &+[ Random(i,R,7) : i in [0..d]]>; 

//some polynomials 
F := randomAf fine (2) ; 
G := randomAf fine (6) ; 
H := randomAf fine (7) ; 

//generators of 1(V(F) \cup V(H,G)) 
I := Ideal([F*G,F*H]); 

/ / experiment 

null := [0 : i in [1. .#Basis(I)]] ; 
t := {**}; 

for j in [1. .700] do 

point := Random (V); 

Include(~t, null eq [Evaluate (Basis (I) [i] .point) : i in [1 . .#Basis(I)] ] ) ; 
end for; 

//result 
t; 

Experiment B.2.4. Count points and tangent spaces on a reducible variety. 
K := FiniteField(7) ; //charakteristik 7 

R<xl ,x2,x3,x4,x5,x6> := PolynomialRing(K,6) ; //6 variables 

V := CartesianPower(K,6) ; 

//reindom affine polynomial of degree d 

randomAffine := func< d I &+[ Random(i ,R,7) : i in [0..d]]>; 

//some polynomials 
F := randomAf fine (2) ; 
G := randomAf fine (6) ; 
H := randomAf fine (7) ; 

//generators of I(V(F) \cup V(H,G)) 
I := Ideal([F*G,F*H]) ; 

null := [0 : i in [1. .#Basis(I)]] ; 
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//the Jacobi-Matrix 

J := JacobianMatrix(Basis(I) ) ; 

size := [NumberOfRows(J) .NumberOf Columns(J)] ; 

A := RMatrixSpace(R,size [1] , size [2] ) ; 
B := KMatrixSpace(K,size[l] ,size[2] ) ; 

t := {**}; 

time 

for j in [1. .700] do 
point := Random (V) ; 

substitude := map< A -> B | x :-> [Evaluate (t, point ) : t in ElementToSequence(x)]>; 
if null eq [Evaluate (Basis (I) [i] , point) : i in [1 . .#Basis(I)] ] 

then Include(~t, Rank(substitude(J) ) ) ; 
else 

Include (~t,-l) ; 
end if ; 
end for; 

//result 
t; 

Experiment B.2.9. 

K := FiniteField(7) ; //charakteristik 7 

R<xl,x2,x3,x4,x5,x6> := PolynomialRing(K,6) ; //6 variables 
V := CartesicOiPower (K,6) ; 

//consider an 5 x 5 matrix with degree 2 entries 
r := 5; 
d := 2; 

Mat := MatrixAlgebra(R,r) ; 



//random matrix 

M := Mat! [Random (d,R, 7) : i in [l..r~2]]; 

//calculate determinant and derivative w.r.t xl 
time F := Determincint (M) ; 
time Fl := Derivative(F, 1) ; 

//substitute a random point 
point := Random (V) ; 
time Evaluate (Fl, point) ; 

//calculate derivative with epsilon 
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Ke<e> := Af f ineAlgebra<K,e I e"2>; //a ring with e*2 = 

Mate := MatrixAlgebra(Ke ,r) ; 

//the first unit vector 
el := <>; 

for i in [1 . .6] do 
if i eq 1 then 

Append(~el ,e) ; 
else 

Append(~el,0) ; 
end if ; 
end for; 

//point with direction 

point 1 := < point [i]+el [i] : i in [1..6]>; 

time Mate ! [Evaluate (x,pointl) : x in ElementToSequence (M)] ; 

time Determinant(Mate! [Evaluate(x,pointl) : x in ElementToSequence (M) ]) ; 

//determinant at 5000 raindom points 

time 

for i in [1. .5000] do 

point := Random (V) ; //random point 

pointl := <point [i] +el [i] : i in [1..6]>; //tangent direction 
//calculate derivative 

_ :=Determinant (Mate! [Evaluate (x, point 1) : x in ElementToSequence (M) ]) ; 
end for; 

Experiment B.4.3. 

R<x,y> := PolynomialRing(IntegerRing() ,2) ; //two variables 
//the equations 

F := -8*x~2-x*y-7*y~2+5238*x-11582*y-7696; 
G := 4*x*y-10*y"2-2313*x-16372*y-6462; 
I := Ideal([F,G]); 

//now lets find the points over F_p 
fimction allPoints (I ,p) 
M := []; 

K := FiniteField(p) ; 

A := Af f ineSpace(K,2) ; 

R := CoordinateRing(A) ; 

for pt in CartesianPower (K,2) do 
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Ipt := IdeaK [R I Evaluate (Basis (I) [k] ,pt) : k in [1 . .#Basis (I)] ] ) ; 

SIpt := Scheme (A, Ipt) ; 

if Codimension(SIpt) eq then; 

Append (~M,pt) ; 
end if; 
end for; 
return M; 
end function; 

for p in PrimesUpTo(23) do 

print p, allPointsd.p) ; 
end for; 

/♦Chinese remaindering 

given solutions mod m and n find 

a solution mod m*n 

soli = [n, solution] 

sol2 = [m, solution] */ 

function chinesePair (soil , sol2) 

n := soll[l] ; 

an := soli [2] ; 

m := sol2[l] ; 

am := sol2[2] ; 

d,r,s := Xgcd(n,m) ; 

//returns d,r,s so that a*r + b*s is 

//the greatest common divisor d of a and b. 

amn := s*m*an+r*n*cmi; 

amn := amn - (Round (amn/ (m*n) ))* (m*n) ; 
if d eq 1 then 

return [m*n,amn] ; 

else 

print "m and n not coprime"; 
return false; 
end if ; 
end function; 

/*take a list {(n_l,s_l) , . . . , (n_k,s_k)} 
and return (n,a) such that 
n=n_l* ... * n_k and 
s_i = a mod n_i*/ 
function chineseList (L) 

//#L >= 2 

erg := L[l] ; 

for i in [2. .#L] do 

erg := chinesePair (L [i] , erg) ; 

end for; 

return erg; 
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end function; 



//x coordinate 

chineseList ( [ [2 , 0] , [5 , 4] , [17 , 10] , [23 , 15] ] ) ; 
//y coordinate 

ChineseList ( [ [2 , 0] , [5 , 1] , [17 , 8] , [23 , 8] ] ) ; 

//test the solution 
Evaluate (F , [1234 , -774] ) ; 
Evaluate (G , [1234 , -774] ) ; 

Experiment B.4.6. Rational recovery, as suggested in von zur Gathen in [23l 
Section 5.10]. Uses the functions allPoints and chineseList from Experiment 
B.4.3. 

R<x,y> := PolynomialRingClntegerRingO ,2) ; //two variables 
//equations 

F := 176*x~2+148*x*y+301*y"2-742*x+896*y+768; 
G := -25*x*y+430*y'2+33*x+1373*y+645; 
I := Ideal ( [F,G] ) ; 

for p in PrimesUpTo (41) do 

print p, allPoints (I ,p) ; 
end for; 

// X coordinate 

ChineseList ( [ [2 , 1] , [13 , 5] , [31 , 7] , [37 , 14] , [41 , 0] ] ) ; 
// y coordinate 

ChineseList ( [ [2 , 0] , [13 , 10] , [31 , 22] , [37 , 18] , [41 , 23] ] ) ; 

//test the solution 

Evaluate (F, [138949,-526048] ) ; 

Evaluate (G, [138949,-526048] ) ; 

/*take (a,n) and calculate a solution to 

r = as mod n 

such that r , s < sqrt (n) . 

return (r/s)*/ 

function recoverQQ (a,n) 



rO 


=a 


sO 


=1 


to 


=0 


rl 


=n 


si 


=0 


tl 


=1 


r2 


=0 



54 



s2:=0; 
t2:=0; 

k := Rouiid(Sqrt(rl*1.0)) ; 
while k le rl do 

q := rO div rl; 

r2 := rO-q*rl; 

s2 := sO-q*sl; 

t2 := tO-q*tl; 

rO : =rl ; 

sO:=sl; 

tO:=tl; 

rl:=r2; 

sl:=s2; 

tl:=t2; 
end while; 
return (r2/ s2) ; 
end function; 

//x coordinate 

recoverQQ (138949, 2*13*31*37*41) ; 

//y coordinate 

recoverQQ (-526048 , 2* 13*31*37*41) ; 

//test the solution 
Evaluate (F , [123/22 , -77/43] ) ; 
Evaluate (G, [123/22,-77/43] ) ; 

Experiment B.4.10. Lifting solutions using p-adic Newtoniteration (as suggested 
by N.Elkies). Uses the function allPoints ftom Example B.4.3. 

//calculate reduction of a matrix M mod n 
function modn(M,n) 

return Matrix (Mrows (M) ,Mcols (M) , [x - Round(x/n)*n : x in Eltseq(M)]); 
end function; 

//divide a matrix of integer by an integer 

//(in our application this division will not have a remainder) 
function divn(M,n) 

return Matrix (Mrows (M) ,Mcols (M) , [x div n : x in Eltseq(M)]); 
end function; 

/ / invert number mod n 
function invn(i,n) 

a,b := Xgcd(i,n) ; 

if a eq 1 then 
return b; 

else return false; 
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end if; 
end function; 

//invert a matrix mod n 
//Ma square matrix over ZZ 

//(if M is not invertible mod n, then is returned) 
function invMatn(M,n) 
Mn := modn(M,n) ; 

MQQ := MatrixAlgebra(RationalField() ,Nrows(M)) !Mn; 
detM : = Determinant (Mn) ; 
if Type ( invn (detM, n)) eq BoolElt then 
return 0; 

else 

return 

(MatrixAlgebradntegerRingO ,Nrows(M)) ! 
(modn(invn(detM,n)*detM*MQQ~(-l) ,n))) ; 

end if; 
end function; 

//(P.eps) an approximation mod eps (contains integers) 
//M affine polynomials (over ZZ) 

//J Jacobian matrix (over ZZ) 

//returns an approximation (P,eps"2) 
function newtonStep(Peps,M, J) 

P := Peps[l] ; 

eps := Peps [2] ; 

JatP:=Matrix(Ncols(J) ,Nrows(J) , [Evaluate(x,Eltseq(P)) : x in Eltseq(J)] ) ; 
JPinv := invMatn(JatP,eps) ; 

MatP:= Matrix (1,#M, [Evaluate (x,Eltseq(P)) : x in Eltseq(M)]); 
correction := eps*modn (divn(MatP*Transpose (JPinv) , eps) , eps) ; 
return <modn(P-correction,eps~2) ,eps~2>; 
end function; 

//returns an approximation mod Peps~(2~num) 
function newton(Peps,M, J,num) 

localPeps := Peps; 

for i in [l..num] do 
localPeps := newtonStepdocalPeps ,M, J) ; 
print localPeps; 

end for; 

return localPeps; 
end function; 

//c.f. example 4.3 

R<x,y> := PolynomialRingdntegerRingO ,2) ; //two variables 
//the equations 
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F := -8*x"2-x*y-7*y~2+5238*x-11582*y-7696; 

G := 4*x*y-10*y~2-2313*x-16372*y-6462; 

I := Ideal([F,G]); 

J := JacobiaiiMatrix(Basis(I) ) ; 

Ap := allPoints(I,7) ; 
for X in Ap do 

MatP := Matrix(Nrows(J) ,Ncols(J) , [EvaluateCj ,x) : j in Eltseq(J)] ) ; 
print x, (0 ne Determinant (MatP) ) ; 
end for; 

Peps := <Matrix(l,2, [2,3] ) ,7>; 
newton(Peps,Basis(I) , J,4) ; 

Peps := <Matrix(l,2, [5,5] ) ,7>; 
newton(Peps,Basis(I) , J,4) ; 

Experiment B.5.1. Count singularities of random surfaces over F5. 

K := FiniteField(5) ; //work in char 5 

A := Aff ineSpace(K,4) ; 

R<x,y,z,w> :=CoordinateRing(A) ; //coordinate ring of IP~3 
//Calculate milnor number 

//(For nonisolated singularities eind smooth surfaces is returned) 
function mu(f ) 

SJ := Scheme (A, (Ideal ( [f] )+JacobianIdeal(f ))) ; 

if 3 eq Codimension(SJ) then 

return Degree (Pro jectiveClosure(S J)) ; 

else 

return 0; 
end if ; 
end function; 

//look at 100 random surfaces 

M := {**>; 

time 

for i in [1 . . 100] do 

f := Random(7,CoordinateRing(A) ,5) ; 

Include(~M,mu(f ) ) ; 
end for; 

print "M:", M; 

Experiment B.5.2. Count singularities of mirror symmetic random surfaces. Uses 
the function mu from Example B.5.1. 



57 



K := FiniteField(5) ; 
A := Af f iiieSpace(K,4) ; 
R<x,y,z,w> :=CoordinateRing(A) ; 



/ /woryi in char 5 
//coordinate ring of IP~3 



//make a random f mirror sjraimetric 
function mysjrmCf) 

return (f + Evaluate (f ,x,-x)) ; 
end function; 

//look at 100 random surfaces 
M := -[**}; 

time 

for i in [1. .100] do 

f :=R!mysym(Raiidom(7,CoordinateRing(A) ,5)) ; 
Include (~M,mu(f) ) ; 
end for; 

print "M:", M; 

Experiment B.5.3. Count ^1-singularities of invariant surfaces. Uses the func- 
tion mu from Experiment B.5.1. 

K := FiniteField(5) ; //work in char 5 

A := Af f ineSpace(K,4) ; 

R<X,Y,Z,W> :=CoordinateRing(A) ; //coordinate ring of IP~3 
RS<x,y,z,w,al,a2,a3,a4,a5,a6,a7> := PolynomialRing(K, 11) ; 

//the 7-gon 

P := X*(X-6-3*7*X"4*Y"2+5*7*X-2*Y-4-7*Y"6) 

+7*Z* ( (X"2+Y"2) -3-2~3*Z-2* (X~2+Y"2) -2+2-4*Z-4* (X~2+Y-2) ) -2-6*Z-7 ; 

//parametrising invariant 7 tics with a double cubic 

U : = (z+a5*w) * (al*z"3+a2*z"2*w+a3*z*w~2+a4*w"3+ (a6*z+a7*w) * (x~2+y"2) ) "2 ; 

//reindom invariant 7-tic 
fimction randomlnvO 

return (P - Evaluate (U, [X,Y,Z,W] cat [Random(K):i in [1..7]])); 
end function; 

//test with 100 examples 
Ml := {**>; 

time 

for i in [1. .100] do 

Include ("Ml, muCrandomlnvO)) ; 

end for; 

print "Ml:", Ml; 
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//singularities of f 
function numAl(f) 

singf := Ideal ( [f] )+JacobiaiiIdeal(f ) ; 
Ssingf := Scheme(A,singf ) ; 
if 3 eq Codimension (Ssingf ) then 
T := Scheme(A,Ideal( [f] )) ; 
//calculate Hessian 
Hess := HessiaiiMatrix(T) ; 
ssf := singf + Ideal (Minors (Hess, 3) ) ; 
if 4 eq Codimension (Scheme (A, ssf ) ) then 

return Degree (ProjectiveClosure (Ssingf)) ; 

else 

return 0; 
end if ; 

else 

return 0; 

end if ; 
end function; 

//test with 100 examples 

M2 := {**}; 

time 

for i in [1 . . 100] do 
Include(~M2,numAl(randomInv())) ; 
end for; 

print "M2:", M2; 

Experiment B.5.4. Estimate number of Al-singularities by looking at y = 0. Uses 
the function numAl from Experiment B.5.3. 

K := FiniteField(5) ; //work in char 5 
A := Aff ineSpace(K,4) ; 

R<X,Y,Z,W> := CoordinateRing(A) ; //coordinate ring of IP~3 
RS<x,y,z,w,al,a2,a3,a4,a5,a6,a7> := PolynomialRing(K, 11) ; 

//the 7-gon 

P := X*(X~6-3*7*X"4*Y"2+5*7*X~2*Y~4-7*Y~6) 

+7*Z* ( (X"2+Y"2) -3-2~3*Z-2* (X~2+Y"2) '2+2~4:*Z~4* (.X~2+Y~2) ) -2-6+Z-7 ; 

//parametrising invariant 7 tics with a double cubic 

U : = (z+a5*w) * (al*z"3+a2*z"2*w+a3*z*w~2+a4*w"3+ (a6*z+a7*w) * (x"2+y"2) ) "2 ; 

//estimate number of nodes 
function numberof singO 
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r := [Random (K) : i in [1..7]]; 

f := Evaluate(P - Evaluate(U, [X,Y,Z,W] cat r),Y,0); 
singf := Ideal ([f])+ JacobianldeaKf ) ; 
Ssingf := Scheme (A, singf ) ; 
if 2 eq CodimensionCSsingf ) then 
//calculate Hessian 
S := Scheme(A,Ideal([f])) ; 
Hess := HessianMatrix(S) ; 
ssf := singf + Ideal (Minors (Hess, 2) ) ; 
Sssf := Scheme(A,ssf ) ; 
if 3 eq Codimension(Sssf ) then 

d := Degree (Pro jectiveClosure (Ssingf )) ; 

//points on the line x=0 

singf X := singf + Ideal ([X]); 

Ssingf X := Scheme (A , singf x) ; 

dx := Degree (Pro jectiveClosure (Ssingf x) ) ; 

if 2 ne Codimension(Ssingf x) then 

dx := 0; 
end if ; 

d3 := (d-dx)*7+dx; 

if d3 ge 93 then 

return <d, d-dx, dx, d3>, r; 

else 

return <d, d-dx, dx, d3>, _; 
end if ; 

else 

return <-l,0,0,0>,_; 
end if ; 
else return <0,0,0,0>,_; 
end if ; 
end function; 

Ml := {**}; 
Mlhit := {**}; 
time 

for i in [1 . . 10000] do 

a,b := number of sing ; 
if assigned (b) then 

Include ("Mlhit, <a,b>) ; 

else 

Include ("Ml , a) ; 

end if ; 
end for; 

print "Ml:", Ml; 
print "Mlhit:", Mlhit; 
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//test 

f := P-Evaluate(U, [X,Y,Z,W,1,2,2,1,1,0,1]) ; 
numAKf ) ; 
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